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I. INTRODUCTION 

The creation of many particle entangled states in 
macroscopic systems is one of the major goals in the stud- 
ies on fundamental aspects of quantum theory |], [|, ||, 
|], H . The notion of entanglement in macroscopic ensem- 
bles allows to investigate the boundary between quan- 
tum physics and classical physics and, possibly, could also 
give some insight into the measurement process [|[ [f[ ||. 
The experimental creation of many particle entangle- 
ment could also lead to the realization of several of the 
"Gedankenexperiments" proposed in the early days of 
quantum theory ||. Also, apart from the fundamental 
physical interest in entanglement, the whole field of quan- 
tum computing and quantum information is based upon 
the ability to create and control entangled states [|l0[ . 

The experimental achievement of atomic Bose-Einstein 
condensation (BEC) Q has opened fascinating possibil- 
ities for studying quantum properties of a macroscopic 
number of cold quantum degenerate atoms in the labora- 
tory. Interesting aspects of many particle entanglement 
can be studied by us ing condensates with internal degrees 
of freedom g| |l| ^ [l| 0, 0, @. For instance, it 
has been shown that the coherent collisional interactions 
in a BEC allow to generate substantial many-particle 
entanglement in the spin degrees of freedom of a two- 
component condensate Q during the free evolution of 
the condensates. 

The extremely long coherence time in a BEC is one 
of the key features making the proposals for engineering 
many particle entanglement in a BEC feasible [[L9| . How- 
ever, one still has to make sure that (i) the creation of the 
many particle entangled state takes place on a time scale 
much shorter than the coherence time, and (ii) the pro- 
duced many particle entangled states are robust against 
decoherence, i.e., particle loss should not destroy the en- 
tanglement. 

In this paper we will investigate in detail a scheme that 
allows to create many particle entangled states in a two- 
component BEC interacting with a classical laser field or 
microwave field on a time scale proportional to 1/N with 
N the number of condensate particles (see also [p0|). Let 
us briefly explain the basic idea of this scheme. If we 



concentrate on the dynamics of the internal states of the 
condensate (e.g. atomic hyperfine levels) the Hamilto- 
nian H of the system is given by (cf. Sec. II A 3| ) 



H = x s 2 z + ns x , 



(i) 



where x 1S determined by the interaction strengths be- 
tween the condensate particles and Q is the Rabi fre- 
quency of the external field interacting with the con- 
densate. Here the angular momentum operator S = 
{S x , S y , S z } is defined by 



S z 



i(a f 6 + 6 f a), 
-(a^a — 6^6), 



where a (b) are bosonic destruction operators for parti- 
cles in internal state A (B) with a fixed spatial mode 
function ipA(B)- We consider the situation, where ini- 
tially all the condensate particles are in internal state \A). 
This state corresponds to an eigenstate of S z with eigen- 
value N/2. Then a tt/2 pulse is applied to the conden- 
sate which brings each particle in a superposition state 
(\A) + \B))/y/2 corresponding to an eigenstate of S x with 
eigenvalue N/2. The internal wave function of the 
condensate is thus given by 



|*(f = 0)> 
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(a f + 6 t ) iV \vac) , 



V2 N N\ 



(2) 



where |vac) is the vacuum state. 

To get a qualitative understanding of the subsequent 
time evolution of \Sfr(t)) we use the familiar phase model 
|l|, H where one replaces S z — > —id§ and S x — ► 
7V/2cos(/>, with <p the relative phase between the two 
condens ate c omponents, in the Hamiltonian and finds 
(cf. Sec. pBl) 



qn 



■ COS ( 



(3) 



We note that the eigenstates of (j) given by 



-imcf) 



\N/2 — i7i)a\N/2 + m)B are entangled states 
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of particles in states A and B |2^, ^3], 24 . Projec- 
tion on these eigenstates gives the phase wave function 
^f(t,(p) = ((p\^f(i)). For large condensate particle num- 
ber N the phase wave function corresponding to the 
initial state \£(0, 0) is well approximated by a narrow 
Gaussian wave-packet with a width of a = 1/y/N cen- 
tered at (j) = 0, i.e. a maximum of the phase potential 
V((/>) = QS cos as schematically shown in Fig. 0a. The 
time evolution due to will first squeeze this wave 
packet due to the harmonic terms of V((j)) resulting in 
an increased width. Then the anharmonicity of V(4>) be- 
comes important and the wave packet splits up as shown 
in Fig. |l|b. These two wave packets correspond to wave 
packets of opposite relative particle number (measured 
by S z ) since they move in opposite direction. Their su- 
perposition is thus a macroscopically entangled state. In 
the remainder of this paper we will investigate the pro- 
cess described above quantitatively. In particular, we will 
use a semiclassical approximation which allows to analyt- 
ically estimate properties of the many particle entangled 
states like the positions of maxima in the relative particle 
number distribution. 

The paper is organized as follows: In Sec. [□] we will 
introduce the model. We derive a two mode Hamiltonian 
describing a two component condensate interacting with 
a laser or microwave field. Then we discuss the phase 
model already used in the introduction and present a 



semiclassical model. In Sec. Ill we define the most com- 



mon types of many particle entangled states and give a 
short overview of recent proposals for entanglement cre- 
ation in BEC's. We present our results in Sec. |^ and 
Sec.0 where we compare the properties of the many par- 
ticle entangled states found by numerical solution of the 
exact Schrodinger equation with the semiclassical esti- 
mates. Furthermore in Sec. [Vl| we improve the scheme 
by one-axis-pre-squeezing. Finally, in Sec. VII we dis- 
cuss the stability of the many particle entangled states 
under decoherence and investigate possible measurement 
strategies to identify the entangled states. Also we study 
the influence of imperfections in the external field on our 
scheme. 



II. MODEL 



In this section we present the model used to study the 
two component condensate interacting with an external 
field. We write down the Hamiltonian for this system and 
use the two mode approximation to obtain a simplified 
description of the system in terms of angular momentum 
operators. Then we introduce the phase model and a 
semiclassical model. 




FIG. 1: Schematic time evolution using the phase- model, (a) 
The initial state is a Gaussian wave packet centered at a maxi- 
mum of the phase potential, (b) The time evolution according 
to Hcfy splits the wave packet into a superposition of two wave 
packets moving in opposite direction. 



A. Hamiltonian 

We consider a two component BEC consisting of N 
atoms in two different hyperfine states \A) and \B) 
coupled by a Raman laser or microwave field (internal 
Josephson effect [^2[ p4|). The Hamiltonian of this 
system is given by 



H = H A + H E 



(4) 



where Ha and Hb describe the two component conden- 
sate and H ext the interaction with the external field. In 
second quantization the terms are given by {% = 1) 



j d 3 r*{ 



**,(5) 



Here ^ = ^fe(r) (k = {A, B}) are bosonic field opera- 
tors that annihilate a particle at position r in the hyper- 
fine state \k). The trapping potential for atoms in state 
k is denoted by Vk = Vfc(r) and their mass is rn. The in- 
teraction strengths are given by Uaa,Ubb and Uab for 
collisions between particles in state A, B and interspecies 
collisions, respectively. The effective Rabi frequency £7r 
is assumed to be positive, real, and position independent, 
i.e., we neglect the momentum transfer induced by the 
external field. The detuning of the field from resonance 
is denoted by A. 
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1. Two mode approximation 

We assume that the spatial degrees of freedom can be 
described using one spatial mode function for each com- 
ponent g| |l§ [2C§ 

*A(r) = a^ A (r) and tf B (r) = 6^ B (r), (6) 

where ^(r) are real normalized wave functions and a, (b) 
are bosonic annihilation operators destroying a particle 
in the internal state A, (B). They obey the usual bosonic 
commutation relations [a,a^] = 1, [b,tf] = 1, [a, 6] = 
and [a, fot] = 0. 



2. Angular momentum representation 

We use the angular momentum operators S = 
(S x ,S y ,S z ) defined in Eq. (||) and the eigenstates of S z 
with eigenvalue n: 



\n) z = 



(a)) N,2+n (M) 



N/2-n 



y/{N/2 + ri)\(N/2-n)\ 



I vac) , 



(7) 



The operators S fulfill the standard angular momen- 
tum commutation relations [Si,Sj] = iSk, with k 
cyclic. From the Heisenberg uncertainty principle we find 

(ASf)(AS])>\\(S k }\\ 



states. To be more precise the continuous eigenstates of 
(j) are given by 



1 

Itt 



+N/2 

E < 

n=-N/2 



-incf) 



\n) z - 



(10) 



Using the relation S z \<p) = —id^ neglecting terms 
of order 1/N and for Q <C Nx/2 we can rewrite the 
Hamiltonian # [^l] in the phase representation as given 
in Eq. (||). The Hamiltonian describes a single ficti- 
tious particle moving as a pendulum [^j, [2^]. The first 
term of can easily be identified as the kinetic energy 
of the particle and the second term is a conservative pe- 
riodic potential. For large N ^> 1 the initial phase wave 
function \I/(0, 4>) = {(j) \^f(t = 0)) is well approximated by 
a narrow gaussian wave-packet 



*(O,0) 



1 



VW+h 



+N/2 

E 

n=-N/2 



N 

N/2-n 



1/2 



Jn<f> 
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-1/4 



e 4 CT 



(ii) 



The width of the Gaussian is a = 1/yN. The phase 
model is valid for Vt <C xiV 01- Since we do not want 
to restrict ourselves to this case we will now derive a 
semiclassical model valid for arbitrary Q. 



3. Two mode Hamiltonian 

Using Eqs. (^|J^) the Hamiltonian reduces (up to a con- 
stant) to 



H = 5S z +xSZ + SlS x , 



(8) 



where S = uoa — oob + A + (uaa — ubb)(N — l)/2, and 
X = ( u aa + ubb — 2uab)/2 with . The single particle 
ground state energies uoa and cjb, the coupling u^i and 
the effective Rabi frequency ft are given by 



v 2 +v 

2m 



um = Uki J d 3 r|^/e| 2 |^l 



(9) 



In the following we will consider the situation, where x > 
and S = 0. 



B. Phase model 

In the introduction we used the phase model [2^] 
to qualitatively explain how the time evolution according 
to H can be deployed to create many particle entangled 



C. Semiclassical model 

The many particle entangled state we are interested 
in are superposition states of two wave packets centered 
at two different relative particle numbers. We want to 
use a semiclassical model to estimate the position of the 
maxima of these wave packets. We assume the collective 
spin S to behave like a classical quantity, which is possible 
if the discreteness of the energy levels is negligible i.e. 
for N ^> 1. Furthermore the semiclassical treatment will 
only be valid as long as interference effects are negligible. 
For the initial state we consider it will turn out that such 
interference effects become important for times t > 2t c , 
where t c defined in Eq. ( J30| ) is the time it takes to create 
the many particle entangled state. 

Under these conditions we replace the spin operator S 
by c-numbers 



N 



(sin 9 cos 0, sin 9 sin 0, cos 9) . 



(12) 



This implies the factorization of expectation values of 
products of operators like e.g. ({S x , S y }} by 2 (S x ) (S y ). 
From the Heisenberg equations of motion for the operator 
S given by 




-x{Sz,S y } 

x{s z ,s x }-ns z 



(13) 
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we obtain for the time evolution of the angles 6, < 



A. Many particle states 



d 
dt 



X N f -ujs'mcj) 
2 \ 2 cos — uj cot 6 cos ( 



(14) 



with uj = 2Vt/xN. The corresponding vector fields are 
shown in Fig. |^ for different values of uj. 

From the continuous Wigner function Wt [^9| defined 
as the Fourier transform of the quantum characteristic 
function W t given by 

WtM) = (2tt)- 2 j V j ' dn'W t (n'^')e +l ^' n - n ^\ 
W t (n'^') = (V(t)\e< n '^' s *)\y(t)), (15) 
we get the initial Liouville distribution on the sphere as 
1 



sine 



:Wi=o(n = 5cosfl,0). (16) 



For the initial state \^f(t = 0)) the Liouville distribu- 
tion is well approximated by a Gaussian of narrow radial 
width a 1/y/N: 



Po(M) =ATe-^ (s - So)2 



2na 2 



e 2 CT 



^2 (cos 2 <9+0 2 ) 



(17) 



We introduce different kinds of many particle states 
which will be used to characterize the states encountered 
in the numerical calculations performed in the subse- 
quent section [v|. 



1. Coherent spin states (CSS) 

The CSS H^] are eigenstates of the angular momentum 
operator 5 ni = ni • S with eigenvalue N/2, where ni 
is a unit vector pointing in the direction (#,</>). These 
states are completely uncorrelated and can be written as 
a separable product of single-particle states given by 



-iS z <f> 



Cm 



e -is v e \+N/2) z 

cos-e-^V+sin-e+^&t 
2 2 



N 



vac) 



cos -e-^' 2 \A) 



sm-e 



\B) 



(20) 



CSS are minimum uncertainty states with (^AS^ 3 ) = 
(ASn 3 ) = N/4, where ni,n2,n3, are unit vectors or- 
thogonal to each other. 



The semiclassical time evolution of the Liouville distri- 
bution is then given by 

Pt(0,$) = / M'amO' j #'P O (0',0') x 

sin 7 7 (18) 

where 0t{0' , (j)'), 4>t(0' \ 4> f ) are the classical trajectories of 
a fictional particle starting at t = from 0' ,<j)' , i.e. solu- 
tions of Eq. (|l4|). The semiclassical expectation value of 
an operator C(S) is given by 



(C(8)) t 



p7T r+ir 

/ dOsmO / ( 

70 7-7T 



c 



N 

— (sin cos <j), sin 6 sin (j), cos ( 



.(19) 



III. ENGINEERING MANY-PARTICLE 
ENTANGLED STATES 

In this section we first define different kinds of many 
particle entangled states and review several schemes for 
engineering entanglement in two component BECs. We 
then study in detail a scheme which allows the production 
of a many particle entangled state in a two component 
BEC at a time scale proportional to 1/N. 



2. Spin squeezed states (SSS) 

The SSS |30|, |3l) are characterized by a reduced vari- 
ance compared to that of the CSS in one of the spin- 
components S n2 or S n3 whereas the variance in the other 
orthogonal component is correspondingly enhanced. The 
amount of squeezing is determined by the squeezing pa- 
rameter £ given by 



. N(AS ni ) 2 
mm * j,, 

(S n2 ) 2 + (S n3 f 



(21) 



If £ 2 < 1 the state of the atoms is non-separable (i.e. 
entangled) as has been shown in |jl2| . The parameter £ 2 
thus characterizes the atomic entanglement, and we refer 
to states with £ 2 < 1 as "spin squeezed states" . 

We also note that this parameter determines the 
amount of noise-reduction in atomic clocks There- 
fore the robust creation [|l2[ and preservation [|22| of SSS 
might prove useful in enhancing the accuracy of these 
atomic clocks. 



3. Maximally entangled states (MES) 

A second type of macroscopically entangled states can 
be written as 



GHZ), = ±= (\Af N + \Bf N ) 
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71 



-N/2) z ) 



2. Dynamical Schemes 



(|0 = O,0 = O) + |0 = 7r,0 = O))(22) 



They are a generalization for N particles of the well 
known GHZ-state of three particles || ||, Q]. As can 
be seen from Eq. (p2|) these states can be written as ei- 
ther a coherent superposition of N particles being in the 
state A and N being in state B or as a coherent superpo- 
sition of two CSS pointing in opposite directions on the 
Bloch-sphere. 



B. Schemes for producing many particle 
entanglement 

There have been several proposals of how to engineer 
many particle entanglement in BECs. Basically these 
schemes can be divided into two different kinds. The first 
possibility is to engineer the ground state of the system to 
be an entangled state. This can be done by appropriately 
manipulating the interaction between the particles |2^, 
|32[ . The entangled state is then created by cooling the 
system to its ground state. The second possibility is to 
control the dynamics of the system such that an initially 
separable state evolves coherently into an entangled state 



1. Thermodynamical Schemes 

For x > an d ^ = the ground state of the Hamilto- 
nian (|l|) is the number-squeezed ("dual fock") Q state 
|0) 2 , with half of the atoms in \A) and the others in \B). 
As first noted by Cirac et al. (26} for x < and \Q\ < N/2 
the ground state of the system corresponds to an entan- 
gled state and a maximally entangled state of the form 
Eq. (|2^ ) is attained for zero coupling (i.e. The 
major obstacle in cooling to the ground state is the small 
energy gap between the ground and the first excited state 
which scales like x- 



One of the simplest ways [20 to obtain a SSS is to 
start from a one component BEC in A (corresponding 
to a CSS \N/2)_ Z ), to apply a fast tt/2 pulse (ft > x^) 
that rotates the CSS on the Bloch sphere by an angle tt/2 
around the ?/-axis aligning it along the x direction (cor- 
responding to \N/2) X ). The subsequent evolution due to 
the Hamiltonian H (with Q = 0) establishes correlations 
among the particles, creating an SSS with a squeezing pa- 
rameter £ « (3/A/') 2 / 3 /2 for large N |3(J on a time scale 
t ?z 2 x 3 1/6 /xA r2/3 . Taking into account the spatial de- 
grees of freedom and inelastic collisions with background 
particles still allows for considerable squeezing of the SSS 
as was more recently shown in fl2|, |l7[ p7[ . As noted by 
M0lmer et al. Q and Castin [p5| one-axis squeezing 
also provides a perfect GHZ-state at a much later time 
t — tt/2x neglecting, however, the spatial degrees of free- 
dom and decoherence processes. 

Even smaller squeezing parameters £ can be obtained 
by engineering two particle interactions resulting in a 
Hamiltonian of the form H[ nt = x($z ~ By) as recently 
proposed in |l8l and thus implementing two- axis- 
squeezing |3(|. In |18| intermediate molecular states of 
two-atoms are used to create a Hamiltonian of the form 
i^int while in (l4| this is achieved by applying a series 
of laser pulses to the two-component BEC. As discussed 
by Law et al. j36| turning on a small coupling Q <C x^ 
can be used to further improve the squeezing properties. 
Gordon et al. [|o) showed numerically that for Q ~ x^ 
Schrodinger Cat states can be obtained on short time 
scales. 

The main task of the remainder of this paper is to in- 
vestigate in detail dynamical schemes which create many 
particle entangled states of the form 



|Cat(D)) JV = 
|Cat( 7 )) JV = 



1 

i 

71 



D \ 

/ z 



(23) 



7T + 7 



,0 



7T — 7 



,0) ,(24) 



on a time scale proportional to 1/N using the Hamil- 
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tonian H^ec (see also [[20|). Eq. (|23|) is a macroscopic 
superposition of eigenstates of the operator S z charac- 
terized by a "distance" D. The other entangled state 
Eq. (24) is a superposition of two CSS separated by an 
angle 7. Note that for D = and 7 = these states 
are not entangled, whereas for D = N and 7 = tt they 
coincide forming an MES. 



the southern hemisphere appear. In contrast to the weak 
coupling regime they are now centered at <j) — as shown 
in Fig. |4|a2,b2. This distribution corresponds roughly to 
a state of the form Eq. ( p4|) as shown in Fig. Qb4. 



B. Distance 



IV. SEMICLASSICAL RESULTS 

In this section we first identify three different coupling 
regimes, for which the evolution of the initial state has 
a qualitatively different behavior. Then we derive ana- 
lytical approximations for the time scales on which the 
many particle entangled states Eq. ( p3|) are formed, and 
the attainable separation D as a function of Q. We find 
that for Q = states close to MES are created on a 

time scale t c = ln(8N)/xN. 



A. Coupling regimes 

As can be seen from Fig. || the vector field ( |l4|) on the 
Bloch-sphere is qualitatively different for 2Vt/xN = uo < 
1 (weak coupling regime) and uo > 1 (strong coupling 
regime). We will thus investigate these two regimes sep- 
arately and also look at the intermediate case uo = 1 
(critical coupling). 

In the weak coupling regime the trajectories passing 
through the maximum of our initial distribution at 9 = 
7r/2,0 = 0, make a full revolution in <\>. The separatrix 
separates the rotational modes of the pendulum from the 
oscillatory ones. The rotational modes of the pendulum 
rotate continuously either clockwise- or counter-clockwise 
(trajectories in the lower/upper half plane of 6 outside 
the separatrix), whereas the oscillatory modes oscillate 
around (j) = — tt/2 (trajectories inside the separatrix), 
as shown in Fig. ||a. In this regime the time evolution 
first squeezes our initial state along the separatrix. Then 
two elongated peaks aligned along the 0-axis centered at 
(j) — ±7r appear, respectively (see Fig. |3|a2-b2). From the 
reduced density distribution 

Psc(M) = J #P t (n,0), (25) 

P t (n = cos (9,0) = sin 0P t (cos (9,0), (26) 

we find good agreement with the Wigner distribution for 
a state of the form Eq. ( p3| ) as shown in Fig. |^a4-b4. 

In the strong coupling regime the trajectories passing 
through # = 7t/2,</> = 0, do neither perform a full revolu- 
tion in <j) nor separate the rotational from the oscillatory 
modes any more as can be seen from Fig. ^|c. In fact, 
there are no rotational modes for Q > x&- However, 
first the initial state again is squeezed like in the weak 
coupling regime (see Fig. [|al,bl). Then two well sepa- 
rated gaussian peaks - one on the northern and another at 



According to the semiclassical time evolution the 
largest distance D of the many particle entangled state 
is equivalent to the largest separation of the separatrix 
at <p = tt in the weak coupling regime and at <\> = in 
the strong coupling regime. By using the conservation of 
energy E(0, 0) = x(N cos(<9) /2) 2 + QN sin(<9) cos(0) /2 = 
QN/2, we find the shape of the separatrix 



COS (/>[#] 



1 — uo 1 cos 2 



sine 



sin 0[(j)] = ^ cos 0±y(| cos 



■ 00+ 1. 



(27) 



It follows that in the weak coupling regime the distance 
D is given by (cf. Appendix pi) 



D = N^uo~(2- 



(28) 



and in the strong coupling regime the angle between the 
two CSS states 7 is 



7 = 2 arcsin(v // ^(2 — uo)). 



(29) 



In the weak coupling regime the maximum distance D 
increases with uo until uo = 1 where D takes its largest 
possible value D = N (see Fig. 5a). Thus for uo = 1 
we obtain a state close to a MES. For larger couplings 
(uo > 1) a superposition state of the form Eq. (|24| ) is cre- 
ated. The angle 7 decreases with increasing uo until at 
uo = 2 no superposition is obtained according to Eq. (|29|). 
Exact numerical calculations show that for uo > 2 macro- 
scopically entangled states of the latter form are still ob- 
tained but with 7 <C 7r. 



C. Time Scales 

We will now focus on the required time t c to create a 
many particle entangled state by using the semiclassical 
time evolution of the wave packets forming the super- 
position states. It is crucial that t c is short compared 
to decoherence times to successfully create many particle 
entanglement. 

The time t c is approximately equal to the time needed 
to "travel" along the separatrix from the point lying at 
the distance a from </> = 0,# = tt/2 to cos# = D/N, i.e., 
the point where the many particle entangled state forms. 
Combining Eqs. dll]27|, |23) one obtains (cf. Appendix 
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Xtc 



21og V / 2N(2 -oj) + V / 2N(2 - oj) - 1 



N^/uA2 
log [8iV(2 - w)] 



CO) 



(30) 



where the last approximation is valid for 2 N (2 — a;) ^> 1. 
We notice that t c scales as 1/N and reaches its mini- 
mum value for uo ~ 1, i.e., near the critical coupling, see 
Fig. 5b. 



V. NUMERICAL RESULTS 

We calculated numerically the solution of the 
Schrodinger equation for the exact model Eq. (||) with 
a moderate number of particles N up to 10 3 . A com- 
parison of the exact results with the approximations of 
Sec. |^ yields good agreement already for N > 50. Also, 
we calculate the overlap of the cat-states with the 
states defined in Eqs. (p3,p3), i.e., their fidelity 



F D = |iv(Cat(D)|^)| 2 



(31) 



Furthermore we find a partial revival of the initial wave 
function at time t ~ 2t c . We characterize this revival by 
the overlap with the initial CSS state R given by 



R=\(*(0)\*(t))\ 2 = \ x (N/2\*(t))\ 



(32) 



Since R turns out to be larger than 1/2 for 0.1 < uo < 
1.9 it is useful for measurement purposes, especially for 
checkin g the coherence of an MES, as described in detail 
in Sec.( gng) . 




a4) 

P(n,t) 



A. Comparison of Semiclassical with Exact Results 



FIG. 3: Comparison of the evolution of the system for 
oj — 0.14 (weak coupling) and N = 50. The exact dis- 
crete Wigner function Wt(n, 0) is plotted in (al-a3) compared 
with the Lioville distribution Pt{n — cos 6, <j>) in (bl-b3) at 
three characteristic times. 1) t — t c /2 showing a SSS; 2) 
t — t c showing the macroscopic superpositon state; and in 
3) t rsj 2t c showing the revival. a4) exact discrete number 
distribution P(n, t) = \ z {n\^(t)) | 2 ; b4) reduced density dis- 
tribution P sc (n, t). 



We calculate numerically the solution of the Liouville 
equation Eq. ( plf ) and compare it with the exact dis- 
crete Wigner function as defined by Leonhardt |2{| for 
< uo < 2. From Figs. |al-b2, g^l-b2 we see excel- 
lent agreement between the two results for N > 50 and 
t < 2t c . Afterwards interference effects become impor- 
tant and the semiclassical model breaks down as shown 
in Figs. [|a3-b3, |]a3-b3. As can be seen by comparing 
Figs. |a4 -b4 and Figs. ||a4-b4 we obtain a revival of the 
initial state of about 55% [90%] for the exact solution in 
the weak [strong] coupling regime, whereas the overlap 
of the semiclassical solution with the initial state is only 
about 10%. 

Next we proceed to check the validity of Eqs. ( p8|]30| ) 
for D and t c . Numerically we find D and t c by looking for 
the maximum F]j(t) > Fo(t). The numerical values for 
D and t c agree very well with the analytical expressions 
Eqs. for N > 50, see Fig. 5a,b. 



B. Critical Coupling 

We see from Fig. 5 that for the critical coupling uo = 1 
we obtain the largest possible D = iV, a short creation 
time t C: the best fidelity Pd, and the largest value for the 
partial revival R. Also the decrease of Fd with increasing 
number of particles N is slowest at the critical coupling. 
In Fig. ^ we plot the time evolution showing the initial 
squeezing followed by the creation of a state close to MES 
and the partial revival of the initial state with R = 79%. 



VI. PRE-SQUEEZING 

We investigate if better results for Fd are obtainable 
by adjusting the initial variances AS y ,AS z by squeez- 
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FIG. 5: Comparison of the exact with the semiclassical re- 
sults, a) Attainable distance D as a function of cj, exact 
solution for N — 16 (dashed curve), N — 128 (dash-dotted 
curve), N = 1024 (dotted curve), and semiclassical result 
(solid curve), b) Time t c as a function of uj: exact solution 
(dashed curves) compared with the semiclassical result (solid 
curves) for N = 16, 128, 1024. c) The maximum fidelity F D 
as a function of uj. d) Revival R for N — 16 (dashed curve), 
N = 128 (dash-dotted curve), and N = 512 (dotted curve). 
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FIG. 4: Same as Fig. § but for uj = 1.87 (strong coupling). 



ing. We choose one-axis-squeezing, which is easily imple- 
mented by turning off the external field after the initial 
7r/2-pulse for a time r. Afterwards a second pulse with 
angle a is applied before the system is evolved accord- 
ing to Eq. (H). We optimize r and a for obtaining the 
largest Fd and restrict ourselves to the most interesting 
case uo = 1. The results are shown in Fig. ||. A sub- 
stantial improvement in the obtainable fidelity Fjj can 
thus be reached by pre-squeezing, (cf. Fig. ||c) which, 
furthermore, decreases much slower with increasing N 
than without pre-squeezing. 



VII. DISCUSSION 



with N = N - 1 and \N : n) z the eigenstate of S z for N 
particles. 

In our scheme, however, we expect the effects of de- 
coherence to be reduced for two reasons, (i) The many 
particle entangled state is created on a very short time 
scale, for which we expect only few particles to be lost, 
(ii) Losing a single particle does not completely destroy 
the entanglement for D < N since 
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N 



D-l 
2 



L N 
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(36) 



A. Stability under Decoherence 

A major problem of creating a macroscopic entangled 
state is decoherence, affecting every realistic system. For 
instance the entanglement in an MES is destroyed by 
losing a single particle, as is easily shown by calculating 



N 



q|GHZ) ^ 
HGHZ^II 



N:N/2) , (33) 



Recently Savage et al. |37| discussed the effects of de- 
coherence due to the presence of non-condensed atoms. 
For low temperatures the primary effect is to introduce 
phase-damping, which was taken into account by a four 
mode model in 0. They showed that decoherence ef- 
fects are negligible for N w 200, which indicates that 
the present scheme could be a robust implementation for 
creating MES. 




-500 500 -500 500 



FIG. 6: Exact probability distribution P(n,t) for N = 10 3 
particles for (a) the initial state, (b) the state created for uj w 
0.14 corresponding to a distance of D = A/2, (c) the GHZ- 
like state obtained for uj = 1, and (d) the state for u) ~ 1.87 
corresponding to an angle of 7 = 2 arcsin(l/2). 
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FIG. 7: Exact number distribution P(n,t) as a function of n 
and t for uo = 1 and A = 2000. 



B. Measurement 

By measuring S z one obtains P(n, t) as shown in Fig. 
This, however, does not guarantee the state to be macro- 
scopically entangled, since a completely incoherent mix- 
ture 

p = \ \D/2) z <23/2| + \ \-D/2) z (-D/2\ (37) 
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10 1 N 10 3 -200 S y 200 



FIG. 8: a) Optimum values for r, a for different values of A. 
b) Fidelity Fd as a function of r and a for A = 400. c) Com- 
parison of Fd without pre-squeezing (dash-dotted curve) , and 
with pre-squeezing (solid curve) as a function of A, where the 
dots mark the numerically calculated values, d) Initial CSS- 
Q-function Q(0,</>) = |(6>,0|^)| 2 for optimal pre-squeezmg 
for N - 400. 



gives the same result for P(n,t). To distinguish between 
the two cases one has either to do density matrix tomog- 
raphy (endoscopy) H|, [39) or measure the purity of the 
system. 

The latter can be realized by time reversal which can 
experimentally be implemented by using the external 
field to apply a series of pulses as shown in Appendix [a|. 
Replacing H by — H after time t c and measuring S x af- 
ter 2t c the two cases can be distinguished. For the pure 
case one obtains always the same outcome +iV/2, while 
for a mixture multiple outcomes are possible. Also the 
partial revival with R> 1/2 could be used to distinguish 
between the two cases without the difficulty of imple- 
menting —H. 



C. Imperfections in the external field 

Experimentally it is possible to adjust the phase /3 of 
the 7r/2 pulses very precisely, while it is much more diffi- 
cult to exactly fulfill the condition of having a pulse are 
equal to tt/2. Therefore we investigate the influence of 
errors in the pulse- area of the initial 7r/2-pulse only. 

The radius of the Bloch-sphere is AT/2, while the width 
of the initial distribution is \/N /2. Thus we expect that 
the required precision in the pulse area scales like 1 / y/N. 
We confirmed this expectation by numerical simulation 
for small N up to N w 10 3 . 
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VIII. CONCLUSION 

In summary, we have investigated schemes which allow 
the creation of macroscopically entangled states with a 
distance < D < N on time scales proportional to 1/N. 
Within a semiclassical approximation we obtained an- 
alytical estimates for D and the time £ c , which are in 
excellent agreement with the numerical results. We also 
showed that the fidelity of these states can be improved 
significantly by one axis pre-squeezing. We estimated the 
effect of imperfections and decoherence, and believe that 
the presented scheme might be used to produce macro- 
scopically entangled states with present state of the art 
technology. 
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The time evolution operator U at time t = Mr after 
M arbitrary rotations is given by 



A I 



U(t) = U ps (r) M = ]J U ps (r) with 



k=l 



U ps (r) = n i?„ i (7,) t e-^ T/J i?„ i (7 J ) = e~ iH ^ . 
j'=i 

By combining two rotations about around axes ni , n2 
(J = 2) we obtain H eS = x(S^ + S^ 2 )/2, since 

U ps ( T ) = e -ixS 2 ni T/2 e -i X Sl 2 T/2 ^ ^ _ X 

(1 - ixSlj/2) « e-'xK+^W 2 (A4) 

to first order in r. Thus H e R = — H can be implemented 
(up to a constant) by choosing ni = (1, 0, 0),ri2 = (0, 1, 0) 
and making use of the identity 



N N + 2 



(A5) 



APPENDIX A: IMPLEMENTING THE 
HAMILTONIAN -H 



APPENDIX B: SEMICLASSICAL 
CALCULATIONS 



First we investigate the time-evolution of H given in 
Eq. (||) for Q ^> %iV. Then we show how a continuous 
sequence of pulses of the external field can be used to 
generate —H which is needed to implement the measure- 
ment process discussed in Sec. VII B. For simplicity we 
consider the case —H(Q = 0) = — x$z- 

For \Q\ ^> |x|iV, one can neglect x$z m @- A pulse 
of angle 7 and phase (3 is characterized by 



/+00 
1 
-00 



\n(t)\dt = 7 > ft(t) = \ ft(t)\e 
and implements the following time-evolution 

Uq(j) = g- i 7[cos(/3)S' a; -sin(/3)S' y ]^ 



i(3 



(Al) 



Its effect corresponds geometrically to a counterclockwise 
rotation around the unit- vector n = (cos(/3), — sin(/3), 0) 
about an angle 7. The inverse is found by Up (7)^ = 
Up(—j) = Up +7r (j). Since a rotation about 7 around the 
z-axis corresponds to three pulses 



R z {l) =e~^ = U (|) t/_ f (7)^ (|) , 
R z {lf =R z {-i) = U (|) ^(7)^ (|) , 



(A2) 



a rotation around an arbitrary vector n = 
(sin 6 cos 0, sin sin 0, cos 0) about an angle 7 is given by 



#n(7) 



= e -'T nS = U- 



- ) R z {-<t>)U^ + , 



(A3) 



In_this appendix we will derive the approximations of 
First we derive the maximal distance D given in 
for a —> 0, then we calculate the time t c given in 
J) needed to create the macroscopic superposition 
states. For simplicity we restrict J^c!) to the intervals 

Therefore we 




< < tt, -7T < <j> < 7T in Eqs. (|27|, 
have to distinguish for the ± sign in the second relation 
of Eq. (^) for the two regimes of interest, i.e. for 00 < 1 
only the + sign holds, whereas for uo > 1 both signs hold. 

For a^OD corresponds to the maximum separation 
along the separatrix for given 0, i.e. 



D = max N cos { 



(Bl) 



Using Eq. fl27| ) we obtain the maximum value forw<l 
at d> = 7r as 



£>= ATcos6>+[0 = tt] = iV^/l - (sm0+[(j) = tt]) 2 
= iV^l- ^±h^.y = Ny/w(2-u). 

Similarly for 1 < uo < 2 we obtain three values for at 
= given by #0 = and the points of maximum sepa- 
ration 6± for the minus sign in Eq. (p7|). The separation 
7 is given by 

7= 0+ - 01 = tt - 2 arcsin (-f - 
= 7r — 2 arcsin(o; — 1) = 2 arcsin ^^6^(2 — o;)^ (B2) 
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Next we calculate the time t c given in Eq. (|30|). Us- 
ing Eq. ([ll]) and the conservation of energy as stated 
in Sec. IV B we obtain dO/dt along the separatrix for 
< < tt/2,0 < <j> < 7T, as 



dO N . N I 



uj — cos^ 
uo sin 9 



■ (B3) 



Integrating this ordinary differential equation gives 



2 



9{t) 



d6' 

6(0) 

2 /•*(*) 



1 — UJ 1 cos z 
sin# 



-1/2 



/ dz [ 

Jz(0) 



dz 1 - z 1 



Ny/u(2-u) 



log 



z(t) 1 + y / l-z(0) 2 /cj(2- J) 
*(0) 1 + ^/l-z(ty/u(2-uj) 



04) 



where = cos 0(t). For a point on the separatrix at 



a distance <r from # = 7r/2,</> = 0, i.e. the peak of the 
initial distribution, z(0) is obtained from the relation 

a 2 = 4(S-S(0)) 2 /7V 2 = 2(1 -sin (9 cos (^) 

= 2cos 2 0/u = 2z(0) 2 /u. (B5) 

The time needed to travel from z(0) = ^/uj/2cf to = 
= ^/oo(2 — cJ), i.e. the point of maximum separation 
of the two number state [gaussian distributions] forming 
the macroscopic superposition state ( |23|]24| ) is then ob- 
tained using Eq. (B4) as 



xt 



■.log 



1 + y/l - <T 2 /2(2 - to) 
y/^/2(2-Lj) 



2 log 



V2(2 - to) /a 2 + y/2{2 - to) /a 2 



N^lu(2-lo) 
Finally using a = l/y/N we obtain Eq. (p0|). 



{B6) 
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